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ABSTRACT 



We describe a new, faster implicit algorithm for solving the radiation hydrodynamics 
equations in the flux-limited diffusion approximation for smoothed particle hydrodynamics. 
J> ■ This improves on the method elucidated in Whitehouse & Bate by using a Gauss-Seidel it- 

C^ - , erative method rather than iterating over the exchange of energy between pairs of particles. 

The new algorithm is typically many thousands of times faster than the old one, which will 
enable more complex problems to be solved. The new algorithm is tested using the same tests 
• performed by Turner & Stone for ZEUS-2D, and repeated by Whitehouse & Bate. 
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1 INTRODUCTION 

Smoothed particle hydrodynamics (SPH) is a Lagrangian method introduced bv lLucvl il977h and Gingold & Mona ghanl il977h that is fre- 
. quently used for simulating astrophysical fluid problems, for example for star and galaxy formation and supernovae. However, it has also 
1 been used outside of astrophysics, for example for the modelling of tsunami and volcanoes (see Monaghan ll992L for a review). 
. {_( | SPH is normally used to model hydrodynamics, with the inclusion of self-gravity in astrophysical contexts. The smoothing length, which 

■ may be variable both in space and in time, enables SPH to naturally adapt its resolution to reflect the local density distribution. However, in 
^ | many astr ophys i cal sit uations radiation transport is also an important element. Various attempts have been made to include radiati on transport 
intoSPH.lLuc^ lr977l) made the first foray into this field, using the diffusion approximation to examine the fissio n of protostars. | Bro okshaw 
Jl985Lll98fil> also used radiative transfer in the diffusion appr oximation^ mod ellin g stars with entropy gradients. lOxlev & Woolfsonl (2003) 
combined SPH with a Monte-Carlo radiation transport method. Viau. S. 1 1995) and Bastien et al. 1 2004) presented an implicit scheme for the 
diffusion approximation. Finally. Iwhitehouse & Batell2004l) developed an implicit scheme for two-temperature (gas and radiation) radiative 
transfer in the flux-limited diffusion approximation. 

In this paper we describe a significant improvement to the method of Whitehouse & Ba3 120041) . namely the implementation of an 
implicit algorithm which is many thousands of times faster, but still models the same physics. Section|2|begins by sum marising the o rigins 
of the flux-limited diffusion approximation, and continues into a brief overview of the explicit SPH equations derived bv lwhitehouse & Bate! 
120041) . before describing i n detail the new method . Se ction|3]describes the tests w e performed to validate this code. These test results should 
be compared with those o f lTurner & Stonel l200ll) and Whiteho ~& Bate! 120041) . 



2 METHOD 

In a frame co-moving with the fluid, and assuming local thermal equilibrium (LTE), the equations governing the time-evolution of radiation 
hydrodynamics (RHD) are 

^+pV-v = 0, (1) 

»-V p + X+ P , (2) 
Dt c 
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p ~5t (p ) = _V ' F ~ Vv:P + 4nKef)B ~~ CKepE > 

p ~5t (p) = ' v ~ 4 ™ ppB + c<ICep£ • (4) 

4^f£) = _V.P-^F (5) 

c 2 Dt\p) c 

(Mihalas & Mih alaJl984tlTurner & Stonj200l[) . In these equations, D/Dt = d/dt + v • V is the convective derivative. The symbols p, e, v 
and p represent the material mass density, energy density, velocity, and scalar isotropic pressure respectively. The total frequency-integrated 
radiation energy density, momentum density (flux) and pressure tensor are represented by E, F, and P, respectiv ely. 

A detailed explanation of the flux-limited diffusion approximation to the above equations is given bv lTurner & Stonel 1200 lh . Here we 
simply summarise the main points. The assumption of LTE allows the rate of emission of radiation from the matter in equations |3|and|4|to 
be written as the Planck function, B. Equations[2|to[5]have been integrated over frequency, leading to the flux mean total opacity x f , and the 
Planck mean and energy mean absorption opacities, k p and k e . In this paper, the opacities are assumed to be independent of frequency so that 
K r = k e and the subscripts may be omitted. The total opacity, x, is the sum of components due to the absorption k and the scattering a. 

The equations of RHD may be closed by an equation of state, specifying the gas pressure, the addition of constitutive relations for the 
Planck function and opacities, and an assumption about the relationship between the angular moments of the radiation field. 

In this paper, we use an ideal equation of state for the gas pressure p = (y - l)up, where u = e/p is the specific energy of the gas. Thus, 
the temperature of the gas is T g = (y — IJpu/R, = u/c v , where p is the dimensionless mean particle mass, R g is the gas constant and c v is the 
specific heat capacity of the gas. The Planck function is given by B = (cr B /n)T*, where cr B is the Stefan-Boltzmann constant. The radiation 
energy density also has an associated temperature T t from the equation E = 4cr B T*/c. 

For an isotropic radiation field P = |£. The Eddington approximation assumes this relation holds everywhere and implies that, in a 
steady state, equation|5|becomes 

F = _J_ V £. (6) 

3*P 

This expression gives the correct flux in optically thick regions, where xp is large. However in optically thin regions where xp ~ * the flux 
tends to infinity whereas in reality |F| < cE. Flux-lim ited diffusion solves this problem by limiting the flux in optically thin environments to 
always obey this inequality. iLevermore & Pomranin j J 198 lh wrote the radiation flux in the form of Fick's law of diffusion as 

F = -DVE, (7) 

with a diffusion constant given by 

D=* (8) 
XP 

The dimensionless function A(E) is called the flux limiter. The radiation pressure tensor may then be written in terms of the radiation energy 
density as 

P= tE, (9) 
where the components of the Eddington tensor, f, are given by 

f= I(i-/)l+I(3/-l)/w, (10) 

where n = WE/\WE\ is the unit vector in the direction of the radiation energy density gradient and the dimensionless scalar function f(E) is 
called the Eddington factor. The flux limiter and the Eddington factor are related by 

f = A + A 2 R 2 , (11) 

where R is the dimensionless quantity R = \VE\/(xpE). 

Equations ITlto l 1 1 I close the equations of RHD, eliminating the need to so lve equation|5| However, we must still choose an expression 
for the flux limiter, A. In this paper, we choose lLevermore & Pomranin j il98lh 's flux limiter 

MR) = — — — — — (12) 
6 + 3R + R 2 

to allow comparison of our results with those of lTurner & Stonel l200lh and lwhitehouse & Bate! 12004^ . 



2.1 The explicit method 

In lWhitehouse & Bate! JiooH) . we described a method by which equations|2|to|4]can be written in SPH formalism. EquationQdoes not need 
to be solved directly since the density of each particle is calculated using the standard SPH summation over the particle and its neighbours. 



SPH with flux-limited diffusion 3 



We define the specific energy of the gas to be u = e/p, and that of the radiation to be f = E/p. The explicit equations in one dimension are 
then 
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(13) 



(14) 



(15) 



where a = 4<tb/c, m, is the mass of SPH particle i, ry = r, - is the difference in positions between particles i and 7, Vy = v, - v^, and 
Wy = W(rij,hij), where W is the standard cubic spline ker nel and the mea n smoothing length of particles i and j is h is = (hj + hj)/2. The 
smoothing lengths are defined in the same manner as those in Whit ehouse & Batj f2004). so each particle has approximately eight neighbours 
unless otherwise stated. 

We use the standard SPH artificial viscosity 



n„ 



(-a v c s ntj +/3 V (it/) /pij 



if v y - r,j > . 



where p,j = h (v, - ■ (r, - / f|r; - + if j, with rf = O.Olh 2 to prevent numerical divergences if particles get too close together. We 
use cr v = 1 and /? v = 2 unless stated otherwise. In equation fT4l the first term on the right hand side is the diffusion term, the second is the work 
done on the radiation field (in one dimension), and the final term allows energy transfer between the radiation and the gas. In eQuation fT5l the 
energy transfer term occurs with opposite sign, while the remaining term is the symmetric SPH expression for work and viscous dissipation 
done on the gas when the thermodynamic variable of integration is energy. 

We also tested a supercritical shock (see section l3~4l with an implicit method derived from the asymmetric variant of the work term in 
equation ll5l 



Diij 
~Dt 



7=1 \Pi ' 



+ acKj 



Pjk 
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(16) 



The two forms of the gas energy equation gave results that differed by two per cent or less. All results presented in this paper use the 
symmetric version. 



2.2 Implicit method 

The implicit method for solving the energy equations described bv lwhitehouse & BatdEool calculated the gas work and viscous terms and 
the diffusion term as an interaction between pairs of particles, subtracting energy from one particle and adding it to another. The radiation 
pressure term was added to £ , while the interaction term between the gas and the radiation was calculated by the solution of a quartic equation. 
The required timestep df was split into N substeps and the particles swept over. This solution was then compared to that with 2N substeps. If 
the fractional error between the two solutions was not less than a specified tolerance, the number of substeps was doubled until the required 
tolerance was reached. 

The new formulation uses a Gauss-Seidel method to ite rate towar ds the solut ion of the system of equations. We use a backwards Euler 
implicit method rather than the trapezoidal method used in Iwhitehouse & Bate] 12004I) . because the former allows larger timesteps to be 
taken. To advance a time-dependent variable A, from time t = n to t = n + 1, the backwards Euler scheme states 

Ar i =A » + di |^ij" +1 . (17) 

For a Gauss-Seidel method involving interactions between particles i and j, the new value A" +1 can be solved for by arranging the implicit 
equations into the form 

A" -dtY i cr i j(A" +l ) 
1 -dfZ/0">j 

where a-,j contains quantities other than A, and A" +i begins as A" and is updated as soon as new values become available. This equation is 
iterated over until convergence is achieved. 

The backwards Euler form of eQuation fT4l is given by 
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C x = % +dtY ^Lbcip^* 1 -Pjtf) ^ -<fc(V • v),m; +1 - Otoe* 
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(21) 



substituting in (y - l)«p for p using our equation of state. These equations can be rearranged into the form of eauation fT8l to solve for 
and u" +[ and Gauss-Seidel iteration performed with all other independent variables fixed. 

We investigated two different approaches to solving these equations. The first method was to iterate over the Gauss-Seidel form of 
equations 1191 and |21] se parate ly but within the same iterations. This resulted in implicit integration that was many times faster than that 
of Whitehouse & Bate 1 2004) in optically thin regimes, but in optically thick regions the performance of the code was similar to that of 
Whitehouse & Bate 1 2004). The method also failed to converge for large timesteps when the energy transfer term (the last terms in equations 
I19l and l2"fl became large (due to high k or large temperature differences between the gas and the radiation). 

By far the most effective method is to solve equations ll9l and [Til simultaneously for m" +i , and to perform Gauss-Seidel iteration on the 
resulting expression. The resulting value of m" +1 is then substituted into equation [T9l to obtain during the same iteration. To simplify the 
subsequent equations we define the following quantities: 
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(V-v) f / i; 
dtD di - dtR pi . 

Using these new variables we can solve equation l2"Tl for 

= \ - «" - df p nj - dr PdX + 1 + r [< +1 ] 4 ) . 



(22) 



The right hand side of equation |22| then replaces in equation lT9l forming a quartic equation in u" +> . If the quartic equation is cast in the 
form «4X 4 + + + a\x + «o = then the co-efficients are given by: 

a 4 = rdrOr-i) 

«3 = 

a 2 = 

a = pe i+ (x-p-i) (-< - dt P nj ) + dt D nj /3 

Solving this quartic equation yields a value for u" +l , which may then be substituted into 



3 



U>; + dt A v + d?r[ ; /' +1 l 4 ) 

•n+l _ V / 



(23) 



using the quantities defined above (for the analytic solution of a quartic equation, see Appendix A of Iwhitehouse & Batell2004i) . These 
solutions for and u" +> are iterated until they converge. 
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2.3 Prediction of position, density and smoothing length 

We found that the accuracy when taking large implicit timesteps could be improved by predicting forward many of the quantities on the 
right-hand sides of equations ll9l and l2"Tl to time t = n + 1. The quantities x, p, and h can be predicted forwards as 

x" +l = x" + At v" 

p'' +l = pi - At pi (V • vf t 

/,;.'+ 1 = }q + At hi (V • v)? . (24) 
The improvement in accuracy was especially apparent in the case of the supercritical shock (see section l3~4l . 



2.4 Convergence criteria 

We define convergence as being when the values of u" +l and £" +l obtained from the m-th iteration satisfy equations 1191 and [Til to a given 
tolerance (with all occurrences of u" +l and on the right-hand sides of these equations being evaluated from iteration m - 1). Thus, for 
example, we iterate until the fractional errors in f given by 

c M - [?; + * 2, ( P ,f r - Pj ff) Si - at (v . v), Mr 1 -'"- 1 -&ac Ki 

are within a certain tolerance, for which we have typically used 1CT 3 . 

In the event that the method fails to converge, we split the timestep into two halves, and begin the iterations again, using the result of 
the first half timestep as the input to the second. If either fails, we split the timesteps by another factor of two and continue this way until 
the system converges, or it reaches some excessively small fraction of the original timestep, at which stage it is no longer computationally 
efficient to continue the calculation. 




2.5 Timestep criteria 

The integration of the hydrodynamic variables requires that the timesteps obey the Courant condition for the hydrodynamic processes. The 
usual hydrodynamical SPH timestep criteria are 

, _ 

c s + hi |V ■ v\i + 1.2 (a v c s + {3 v hj |V • v|;) 

and 

dfforce,,- = <T a/A' ( 2? ) 

\ kl 

where we use a Courant number of £ = 0.3, unless otherwise noted, and a, is the acceleration of particle i. The lesser of these two quantities 
gives the hydrodynamical timestep. 

There is also an explicit timestep associated with the radiation hydrodynamics, described in detail in lwhitehouse & Batel i20o4) . This 
timestep is typically much smaller than the hydrodynamical timestep, and the new implicit method enables us to forgo the use of smaller 
timesteps in favour of the large hydrodynamical timestep. 



3 TEST CALCULATIONS 

We have once again duplicated the tests done bv lTurner & Stonel i200lh . as we did in lwhiteh ouse & Batel l20 04h. th is time however using 
the new code. This code achieves the same or better accuracy than the code described in Whitehouse & Bate 1 2004). However, in the vast 
majority of cases (especially those involving moving fluids) the new code is significantly faster than the old one. 

3.1 Heating and cooling terms 

We tested the interaction between the radiation and the gas to check that the temperatures of the gas and the radiation equalise at the correct 
rate when T g + T r initially. A gas in a domain 10 cm long with 100 particles was set up so that there was no velocity, with a density p = 10~ 7 
g cnT 3 , opacity k = 0.4 cm 2 g -1 , and y = 4 and E = £p = 10 12 ergs cirT 3 . Two tests were carried out, one where the gas heated until it 
reached the radiation temperature, and one where it cooled. The first test had e = up = 10 2 ergs cnT 3 , and the second e = up = 10 10 ergs 
cirT 3 . The boundaries of the calculation used reflective ghost particles. 

This problem can be approximated in the case where the energy in the radiation is much greater than that in the gas by the differential 
equation 
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Figure 1. The evolution of the gas energy density e as it equilibriates with a radiation energy density E = 10 12 erg cnT 3 . In the upper case, e = 10 10 erg cm" 3 
initially, while in the lower case e = 10 2 erg cm~ 3 . The solid line is the analytic solution, the crosses are the results of the SPH code using implicit timesteps 
of the lesser of 10~ 14 s or five percent of the elapsed time, and the squares with a timestep of the lower of 10~" s or five percent of the elapsed time. The 
symbols are plotted every ten timesteps. 



and assuming E is constant. 

In figure Q the solid line is this analytic solution, plotted both for the cases where T g increases and decreases. The crosses are the 
results of the S PH code us ing an implicit timestep that is set to the greater of 1(T 14 s or five percent of the time elapsed, similar to the way 
Whit ehouse & Bate) 12004) performed this test. The squares are similar, but with a timestep being the greater of 10" 11 s or five percent of the 
time elapsed. As can be seen, the match between the analytic solution and the solutions given by the SPH code is once again excellent. On 
this test, the new method and that of Whitehouse & Bate ( 2004) are of comparable speed and each takes a few minutes. 

3.2 Propagating radiation in optically thin media 

In the standard diffusion approximation, in optically thin regions, radiation propagates at near infinite speed. This is unphysical, and the flux 
limiter has been introduced to limit the diffusion of radiation to the speed of light. To examine how well our code limits the speed of the 
radiation, a one centimetre long one-dimensional box is filled with 100 equally spaced SPH particles, with E = 10~ 2 erg cm 3 = 0.4 erg 
g -1 ), p = 0.025 g cm 4 , and k = 0.4 cm 2 g~' . Initially, the radiation and gas are in thermal equilibrium. 

At the start of the simulation, the radiation energy density for the leftmost ten particles was changed to E = 0.1 erg cm~ 3 (£ = 4 erg 
g -1 ) causing a radiation front that was allowed to propagate across the region. The ghost particles were reflective except in specific radiation 
energy which was fixed equal to £ = 4 erg g -1 at the left hand boundary and £ = 0.4 erg g _I at the right hand boundary. The implicit code 
was used with various timesteps ranging from an explicit timestep to a single implicit step lasting 10~" s. The results are shown in figure|5| 

As shown by the figure, the radiation pulse propagates at the corr ect speed, even using one si ngle large timestep. The front is smoothed 
out in a manner similar to the results of lTurner & Stone! 1200 ll) and [ whitehouse & Bate! 120041) : both methods are quite diffusive in this 
situation. The new code is slower than that of Whitehouse & Bate ( 2004) for an explicit timestep, but superior for all longer timesteps. 

3.3 Optically-thick (adiabatic) and optically-thin (isothermal) shocks 

A shock-tube test, identical to the one in Iwhitehouse & Bate! 12004 ), was set up to investigate the way the code simulated optically-thin 
and optically-thick regimes and the transition between them. In the limit of high optical depth, the gas cannot cool because the radiation is 
trapped within the gas; thus the shock is adiabatic. An optically-thin shock, on the other hand, is able to efficiently radiate away the thermal 
energy and thus behaves as an isothermal shock. In these shock tests, the gas and radiation are highly coupled and, thus, their temperatures 
are equal. 




(28) 
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Figure 2. The propagation of a radiation pulse across a uniform medium. The time is t = 10 s. The vertical dashed line shows the expected position of the 
pulse based on the speed of light. The results are almost independent of the size of the implicit timestep used. Results are given for implicit timesteps equal 
to (solid line), ten times (dotted line), one hundred times (short-dashed line), and one thousand times (long-dashed line) the explicit timestep. The dot-dashed 
line gives the result using a single implicit step of 10~" s. The initial conditions are also shown as a solid line. 



A domain 2 x 10 15 cm long extending from x = -1 x 10 15 to 1 x 10 15 cm was set up, with an initial density of p = 10~ 10 g cm" 3 , and 
the temperatures of the gas and radiation were initially 1500 K. One hundred particles were equally spaced in the domain, with those with 
negative x having a velocity equal to the adiabatic (y = 5/3) sound speed vo = c s = 3.2 x 10 5 cm s _I , and those with positive x travelling at the 
same speed in the opposite direction. The two flows impact at the origin, and a shock forms. Opacities of k = 40, 0.4, 4.0 x 10~ 3 and 4.0 x 10~ 5 
cm 2 g -1 were used to follow the transition from adiabatic to isothermal behaviour. Ghost particles were placed outside the boundaries and 
maintain the initial energies of their respective real particles. The boundaries moved inwards with the same velocity as the initial velocities 
of the two streams. These moving boundaries cause slight perturbations in the densities of those particles closest to the boundaries, however 
this does not affect the solution in the vicinity of the shocks. 

The optically thick and thin limits can be solved analytically (e.g. lZerdovich & RaizeJ2002l) . The shock speed is given by 



(Teff - 3) + yj(y<£ + l) 2 v 2 + 16y c , 



D = , (29) 

where y eff = 1 for the isothermal case, and y cff = 5/3 for the adiabatic case. The ratio of the final to the initial density is given by 

^ = l + £. (30) 
Po D 

and, for the adiabatic shock, the ratio of the final to the initial temperature is given by 

Il = PL + v j£El. (31) 
To Po Po 

These analytic solutions are shown by the solid and dashed lines in figure|3| In the figure, the opacity decreases from top to bottom showing the 
transition from optically-thick (adiabatic) to optically-thin (isothermal) behaviour. The extremes are in good agreement with their respective 
analytic adiabatic and isothermal solutions. Note that the spike in thermal energy near the origin and the corresponding reduction in density 
for the optically-thick case (due to 'wall-heating') is softened by the radiation transport that occurs in the intermediate opacity calculation 
with k = 0.4. 

In comparison to the code described in Whitehouse & BateN2004 these shocks were many times faster. The timestep was limited by 
the hydrodynamical criteria only. The k = 40 and k = 0.4 shocks ran in less than a minute, compared to the old code which took thirty-four 
hours and thirty-five minutes, respectively. The k = 4 x 10~ 3 shock took just under one minute, compared to forty-five minutes for the old 
code. The k = 4 x 10~ 5 code took twenty-three minutes, compared to the previous time of ten hours. 
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Figure 3. A set of shocks with differing opacity at time t = 1.0 X 10 9 s. Density is on the left, and gas temperature on the right. The crosses are the SPH 
results; the solid line gives the analytic solution for an adiabatic shock, and the dashed line for an isothermal shock. The opacities are (top) k = 40, (upper 



middle) k = 0.4, (lower middle) k 
isothermal behaviour. 



4.0 X 10 and (bottom) k = 4.0 X 10 cm g . As the opacity is decreased, the shocks transition from adiabatic to 



3.4 Sub- and super-critical shocks 

A supercritical shock occurs when the photons generated by a shock have sufficient energy to preheat the material upstream. The characteristic 

temperature profile of a supercritical shock is where the temperature on either side of the shock is similar, ra ther than the downstream 

I ■ I 

temperature being much higher than that upstream, as occurs in a subcritical shock (see Zel'dovich & Raizer 2002, for more details). 

The initial conditions of this problem are those of lSincell. Gehmevr & Mihalasl Jl999h and hurner & StoneN200ll) . A gas with opacity 
k = 0.4 cm 2 g" 1 , uniform density p = 7.78 X 10 -10 g citT 3 , mean molecular weight yu = 0.5 and y = | is set up with £ and u in equilibrium, 
with a temperature gradient of T = 10 + (l5x/7 x 10 10 ) K. Initially the particles are equally spaced between x = and x = 1 X 10 10 cm for 
the supercritical shock, and between x = and x = 3.5 x 10 10 cm for the subcritical shock. At time / = a piston starts to move into the fluid 
from the left-hand boundary (simulated by moving the location of the boundary). For the subcritical shock the piston velocity is v p = 6 km 
s~' , and for the supercritical shock v p = 16 km s~' , as per lsincell et alJ il999T> . The ghost particles are reflective in the frame of reference of 
the boundary. Artificial viscosity parameters a v = 2 and /3 V = 4 were used to smooth out oscillations. 

The results of calculations for a sub-critical shock (piston velocity v p = 6 km s~') are shown in figure [4] The top left panel for each 
shows the temperature of the radiation field (solid line) and the gas (dotted line) against position, and the top right shows the same quantities 
against optical depth t, with r = set at the shock front (measured from the density distribut ion). The bottom lef t pa nel shows n ormalised 
flux, and the bottom right the value of the Eddington factor. The analytic solutions discussed bv lsincell et aljfl999l) and lZel'dovich & RaizeJ 
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Figure 4. The sub-critical shock with piston velocity 6 X 10 5 cm s _1 and 100 particles. The large panels show the results using the explicit code. The top 
panels show radiation (solid line) and gas (dotted line) temperatures. The bottom left panel shows the normalised flux and bottom right panel the Eddington 
factor. The long-dashed lines give the analytic solutions for the gas temperature and normalised flux. An Eddington factor of 1/3 is also indicated for reference 
(short-dashed line, lower right panel). The subpanels plot the logarithm of the difference between the results using the explicit code and the implicit code (see 
the main text). The subpanels show the results with (lower subpanel) and without (upper subpanel) predicting x,h and p forwards in time. The implicit code 
was run with timesteps one times (solid line), one tenth of (long-dashed line) and one hundredth of (short-dashed line) the hydrodynamic timestep. 



(2002) for the temperatures and fluxes of the shocks are shown with long-dashed lines. Figures|4]and[5]are plotted using the explicit code. 
In subpanels beneath the main panels, we compare the results from the implicit code with the explicit results. Calculations were performed 
using the implicit code with timesteps of one, one tenth and one hundredth the hydrodynamical timestep criteria. In the subpanels, we plot 
the differences of the implicit results with respect to the explicit results. We divide the difference between the implicit and explicit values 
by the explicit value to obtain a fractional error and take the logarithm of the absolute value of this fraction. Thus, a difference of -2 on the 
subpanels corresponds to an error of 1 percent with respect to the explicit result. The subpanels in figure|4]shows the log fractional error of 
an implicit code with (bottom subpanels) and without (top subpanels) the prediction of x, h and p discussed in section l2~3l Similarly, figure[5] 
show the super-critical shock in the same way, excluding (top) and including (bottom) the prediction mentioned in section l2~3l The prediction 
of x, h and p makes the sub-critical shock significantly more accurate for hydrodynamical timesteps, whilst having a smaller benefit for the 
super-critical shock. Figure|6|shows how increasing the resolution of the simulation enables us to resolve the spike in gas temperature at the 
shock front in a super-critical shock. 

The new implicit method is many times faster than the old method. In lwhitehouse & Bate! |2004) . the super-critical shock with the 
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Figure 5. The super-critical shock, with piston velocity 1.6 X 10 6 cm s -1 , 100 particles. The subpanels shown the results with (bottom) and without (top) 
predicting xji and p forwards in time. This shock is strong enough for radiation from the shock to preheat the gas upstream. See figure|4]for details of the line 
meanings. 



hydrodynamic timestep criteria took nearly ten days. The new code performed the same calculation on a comparable CPU in two minutes, 
making the new algorithm approximately 10 4 times faster. A similar improvement in performance can be seen with the sub-critical shock - 
the hydrodynamic timestep run took over twenty-three days for Whitehouse & Bate 1 2004), while the present method ran in four minutes, 
yielding an increase again of « 10 4 times. 



3.5 Radiation-dominated shock 

In material of high optical depth the radiation generated in a sh ock cannot diffuse away at a high rate, and so the radiation becomes confined 
in a thin region adjacent to the shock. Frurner & Stonel 1200 lb performed a calculation that tests whether the shock thickness is what one 
would expect in these circumstances. An extremely high Mach number shock (Mach number of 658) is set up, with the gas on the left having 
an initial density of p = 0.01 g cm 4 , opacity k = 0.4 cm 2 g~', temperature T r = T g = 10 4 K, and speed 10 9 cm s _1 . The gas on the right has 
density p = 0.0685847 g cm" 3 , opacity k = 0.4 cm 2 g~\ temperature T c = T g = 4.239 x 10 7 K, and speed 1.458 x 10 s cm s~'. The density 
contrast is set initially by using different mass particles, however as the simulation evolves particles from the left enter the shock and by the 
time the figure is plotted all the particles shown have the same mass. The locations of the boundaries move with the same speed as their 
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Figure 6. The super-critical shock, with piston velocity 16 km s , showing how changing the resolution affects the spike in gas temperature at the shock. All 
calculations use a hydrodynamical timestep and the implicit code. The solid line is with 100 particles, the dotted line with 200 particles, and the long-dashed 
with 500 particles. The short-dashed line shows the results with 500 particles and double the number of neighbours (16 instead of 8). 



respective particles, and the properties of the ghost particles outside these boundaries are fixed at their initial values. We use hydrodynamical 
timestep with a Courant number of ( = 0.03. 1500 particles are equally spaced over a domain extending from x = -6x 10 5 cm to x = 1.5x10 s 
cm initially, with the discontinuity at x = 0.5 x 10 5 cm. The location of the shock should be fixed in this frame, although individual particles 
flow through the shock. 

After a period where a transient feature forms at the shock front and drifts downstream with the flow, a stable shock is established. Its 
thickness is expected to be roughly equal to the distance / = cA/Kpui, where u l is the speed of material flowing into the shock front. 

Figure shows the results at a time t = 5 x 10~ 4 s, of the radiation energy density E, gas energy density e, velocity v, and density p 
versus position. The vertical dashed lines indicate the expected shock thickness and the SPH results are in good agree ment. The a fter-effects 
of the transient moving downstream can again be seen on the right of the plots of density and gas energy density. Whitehou s~& Bate! |2004 ) 
were unable to run this test case with the large timestep used here because it required an unacceptable amount of computational time. The 
calculation presented here required approximately a week. 



4 CONCLUSIONS 

We have presented a more efficient method for performing radiative transfer in the flux-limited diffusion approximation within the SPH 
formalism. This gives a speed increase of many thousand times over the code of lwhitehouse & Bate! ^2004) . In every test, the new implicit 
code is much faster than the old code for large implicit timesteps, with no loss in accuracy. 

Whilst the method described here is presented in one dimension, the addition of the algorithm into a three-dimensional code is easily 
accomplished. The major difference between the one- and three-dimensional algorithms is the form of the radiation pressure, which involves 
a more complicated tensor equation. We are performing simulations in three dimensions using this algorithm which will be published in due 
course. 
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Figure 7. The radiation-dominated shock at time ( = 5x 1CT 4 s. We plot the velocity, radiation and gas energy densities and gas density versus position. The 
vertical dashed lines show the expected shock thickness. Gas flows into the shock from the left. The after effects of the transient that occurs at the start of the 
calculation can be seen at the far right of the plots. 



REFERENCES 

Bastien P., Cha S.-H., Viau S., 2004, in Revista Mexicana de Astronomia y Astrofisica Conference Series SPH with radiative transfer: 

method and applications, pp 144-147 
Brookshaw L., 1985, Proceedings of the Astronomical Society of Australia, 6, 207 
Brookshaw L., 1986, Proceedings of the Astronomical Society of Australia, 6, 461 
Gingold R. A., Monaghan J. J„ 1977, MNRAS, 181, 375 
Levermore C. D., Pomraning G. C, 1981, ApJ, 248, 321 
Lucy L. B., 1977, AJ, 82, 1013 

Mihalas D., Mihalas B. W., 1984, Foundations of Radiation Hydrodynamics. Oxford University Press 

Monaghan J. J., 1992, Ann. Rev. Astron. Astrophys., 30, 543 

Oxley S., Woolfson M. M., 2003, MNRAS, 343, 900 

Sincell M. W., Gehmeyr M., Mihalas D., 1999, Shock Waves, 9, 391 

Turner N. J„ Stone J. M., 2001, ApJS, 135, 95 

Viau, S. 1995, Doctoral thesis, University of Montreal 

Whitehouse S. C, Bate M. R., 2004, MNRAS, 353, 1078 

Zel'dovich Y. B., Raizer Y. P., 2002, Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena, edited Hayes and Prob- 
stein. Dover, New York 



